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Fully  three-dimensional  numerical  solutions  can  quantify  exterior  seismic  or 
acoustic  propagation  throughout  complex  geologic  or  atmospheric  domains.  Re¬ 
sults  from  impulsive  sources  typically  reveal  propagating  waves  plus  reverbera¬ 
tions  typical  of  multi-path  scattering  and  wave-guide  behavior,  with  decay  to¬ 
ward  quiescent  motions  as  the  dominant  wave  energy  moves  out  of  the  domain. 

Because  such  computations  are  expensive  and  yield  large  data  sets,  it  is  advan¬ 
tageous  to  make  the  data  reusable  and  reducible  for  both  direct  and  reciprocal 
simulations.  Our  objective  is  efficient  time-domain  simulation  of  the  wave-field 
response  to  sources  with  arbitrary  time  series.  For  this  purpose  we  developed  a 
practical  and  robust  technique  for  superstable  model  identification.  A  supersta¬ 
ble  model  has  the  form  of  a  state-space  model,  but  the  output  matrix  contains  the 
system  dynamics.  It  simulates  propagation  with  the  fidelity  of  the  pulse  response 
calculated  for  the  numerical  system.  Our  development  of  the  superstable  tech¬ 
nique  was  motivated  by  our  initial  application  of  the  Eigensystem  Realization 
Algorithm  to  wave-field  systems  from  high-performance-computing  analyses, 
where  we  recognized  exterior  propagation  features  allowing  superstable  model 
assignment.  Most  importantly  the  pulse  response  and  its  decay  over  the  domain 
are  captured  in  a  finite  duration,  and  decay  to  zero  beyond  a  finite  number  of 
time  steps  implies  a  system  with  zero  eigenvalues.  We  demonstrate  propagation- 
system  identification  with  pulse  response  data  derived  from  supercomputer 
analysis,  and  conclude  that,  using  superstable-identified  systems,  we  are  able  to 
create  reusable  and  reducible  propagation-system  models  that  accurately  simu¬ 
late  the  wave  field  using  a  fraction  of  the  original  computational  resources. 

INTRODUCTION 

Simulating  the  propagation  of  seismic  and  acoustic  energy  involves  solutions  of  wave  equa¬ 
tions.  Fully  three-dimensional  numerical  solutions  can  quantify  propagation  throughout  complex 
geologic  or  atmospheric  domains,  and  are  sometimes  necessary  to  support  a  physical  understand¬ 
ing  of  practical  scenarios  without  compromising  realism.  For  example,  sound  or  vibration  propa- 
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14.  ABSTRACT 

Fully  three-dimensional  numerical  solutions  can  quantify  exterior  seismic  or  acoustic  propagation 
throughout  complex  geologic  or  atmospheric  domains.  Results  from  impulsive  sources  typically  reveal 
propagating  waves  plus  reverberations  typical  of  multi-path  scattering  and  wave-guide  behavior,  with 
decay  toward  quiescent  motions  as  the  dominant  wave  energy  moves  out  of  the  domain.  Because  such 
computations  are  expensive  and  yield  large  data  sets,  it  is  advantageous  to  make  the  data  reusable  and 
reducible  for  both  direct  and  reciprocal  simulations.  Our  objective  is  efficient  time-domain  simulation  of 
the  wave-field  response  to  sources  with  arbitrary  time  series.  For  this  purpose  we  developed  a  practical 
and  robust  technique  for  superstable  model  identification.  A  superstable  model  has  the  form  of  a 
state-space  model,  but  the  output  matrix  contains  the  system  dynamics.  It  simulates  propagation  with  the 
fidelity  of  the  pulse  response  calculated  for  the  numerical  system.  Our  development  of  the  superstable 
technique  was  motivated  by  our  initial  application  of  the  Eigensystem  Realization  Algorithm  to  wave-field 
systems  from  high-performance-computing  analyses  where  we  recognized  exterior  propagation  features 
allowing  superstable  model  assignment.  Most  importantly  the  pulse  response  and  its  decay  over  the  domain 
are  captured  in  a  finite  duration,  and  decay  to  zero  beyond  a  finite  number  of  time  steps  implies  a  system 
with  zero  eigenvalues.  We  demonstrate  propagationsystem  identification  with  pulse  response  data  derived 
from  supercomputer  analysis,  and  conclude  that,  using  superstable-identified  systems,  we  are  able  to  create 
reusable  and  reducible  propagation-system  models  that  accurately  simulate  the  wave  field  using  a  fraction 
of  the  original  computational  resources. 
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gation  from  single  or  multiple  sources,  in  exterior  settings  with  realistic  natural  terrain,  natural 
media,  or  man-made  features,  can  be  visualized  using  results  of  3D  high-performance  computing 
(HPC).  Time-domain  results  typically  reveal  propagating  waves  plus  reverberations  typical  of 
multi-path  scattering  and  wave-guide  behavior,  with  decay  toward  quiescent  motions  after  the 
dominant  wave  energy  from  duration-limited  sources  moves  out  of  the  domain.  Because  such 
computations  are  expensive  and  yield  data  sets  with  a  very  large  number  of  outputs,  i.e.,  beyond 
106  output  locations  in  recent  work,  it  is  advantageous  to  identify  a  reusable  and  reducible  system 
for  both  direct  and  reciprocal  simulations. 


Our  HPC  computations  include  finite-difference  time-domain  (FDTD)  analyses  of  linear 
propagation  media.  For  example,  our  FDTD  computations  for  propagation  within  a  time-invariant 
acoustic  model  solve  the  following  first-order  partial  differential  equations:1 
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where  P(x1,x2,x3,t)  is  pressure,  our  output  variable  of  primary  interest;  v(x,,x2,x3,t)  is  the  parti¬ 
cle  velocity  vector  with  components  v, ,  v2 ,  and  v3;  p(xvx2,x3),  c{xx,x2,x3),  and  cr(xl,x2,x3)  are 
the  spatially  varying  material  density,  speed  of  sound,  and  flow  resistivity,  respectively;  and  Q  is 
a  dilatation-rate  source.  The  method  discretizes  the  solution  domain  using  second-order  finite  dif¬ 
ferences  both  in  time  and  on  a  Cartesian  variable  staggered  grid.  It  applies  an  absorbing  boundary 
condition  at  the  model  edges.2  Leapfrog  computations  solve  the  difference  equations  explicitly, 
starting  from  a  quiet  initial  condition.  We  assign  zero  flow  resistivity  to  the  acoustic  propagation 
media  for  which  we  will  identify  a  linear  system,  and  nonzero  flow  resistivity  to  model  sound 
absorption  by  porous  ground. 

Our  objective  is  efficient  simulation  of  the  wave-field  response  to  sources  with  arbitrary  time 
series.  We  approach  this  as  a  post-processing  step  for  FDTD  analyses,  by  quantifying  the  decay¬ 
ing  pulse  response  over  2D  and  3D  sub  domains  of  our  numerical  linear-time-invariant  (LTI) 
wave-field  system.  Working  with  the  pulse  response  is  a  classical  approach,  but  we  wish  to  make 
use  of  state-space  mathematics  for  its  computational  advantages.  This  requires  identification  of  a 
state-space  model. 

Our  initial  identification  technique  applied  a  variation  of  the  Eigensystem  Realization  Algo¬ 
rithm  (ERA).3'4  We  followed  this  process: 

(1)  By  FDTD  analysis,  generate  responses  over  output  wave-field  domains  using  a  single  in¬ 
put.  The  input  time  series  can  be  impulsive  with  a  dominant  bandwidth  that  extends  no  higher 
than  the  accurate  frequency  range  of  the  FDTD  computation.  Select  the  duration  to  ensure  that 
responses  over  the  output  domains  decay  everywhere  to  an  accepted  quiet  or  still  state,  i.e.,  to 
ensure  an  acceptable  dynamic  range. 

(2)  By  Fourier  analysis,  calculate  frequency  and  pulse  responses  from  input  and  output  time 
series.  Use  the  pulse  response  terms  as  the  Markov  parameters  in  the  identification  problem. 

(3)  By  a  modified-ERA  analysis,5'6  identify  a  state-space  model.  The  model  can  be  used  for 
subsequent  wave-field  simulations  using  inputs  with  the  original  source  location  and  configura¬ 
tion  but  with  user-definable  time  series. 
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Figure  1.  Example  result  using  ERA.  (a)  Geometry  of  model  for  urban  acoustic  propagation  by 
FDTD  analysis.  An  impulsive  sound  source  at  the  model  center,  just  above  ground,  was  used  to  gen¬ 
erate  Markov  parameters,  (b)  Result  from  FDTD  verification  analysis,  showing  snapshot  of  pressure 
field  (Pa)  above  ground  and  rooftops  at  time  =  1.085  s,  from  a  source  containing  early  pulses  and 
continuing  sinusoids.  Triangle  symbol  shows  location  of  extracted  signals,  (c)  Spectragram  and  ex¬ 
tracted  signal  from  FDTD  analysis,  (d)  Corresponding  spectragram  and  signal  from  simulation  using 
ERA-identified  model.  The  ERA-model  simulation  produced  the  illustrated  level  of  fidelity  through¬ 
out  the  pressure  wave  field  using  a  circa  2009  laptop  and  approximately  one  millionth  of  the  re¬ 
sources  of  the  HPC  analysis,  as  quantified  by  (#coresxdurationxmemory  used). 

This  modified  ERA,  similar  to  ERA  using  data  correlations,7  defines  Hankel  matrix  products 
containing  auto-  and  cross-correlation  sequences  of  the  pulse  response.  We  are  specifically  inter¬ 
ested  in  simulation-generated  systems  where  the  number  of  outputs  is  much  larger  than  the  num¬ 
ber  of  inputs,  and  developed  the  modified  ERA  to  require  singular  value  decomposition  (SVD)  of 
a  Hankel  matrix  product  whose  row  and  column  dimensions  depend  on  the  number  of  time  steps 
in  the  pulse  response,  and  not  on  the  number  of  outputs.  For  our  bandwidths  and  durations  the 
number  of  time  steps  is  typically  in  the  hundreds,  and  so  the  SVD  is  of  a  small  matrix.  The  modi¬ 
fied  ERA  successfully  creates  state-space  models  with  model  order  equal  to  the  number  of  time 
steps  of  a  fully  decayed  wave-field  pulse  response,  and  it  supports  subsequent  model-order  reduc¬ 
tion.5'6  Figure  1  highlights  an  example  result. 

From  our  work  with  ERA,  we  recognized  propagation  model  features  allowing  development 
of  what  we  call  a  superstable-system  identification  technique.8  These  features  are  that  the  pulse 
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response  over  the  domain  and  its  decay  to  a  negligible  level  are,  for  practical  purposes,  captured 
in  a  finite  duration;  that  decay  to  zero  beyond  a  finite  number  of  time  steps  implies  a  system  with 
zero  eigenvalues;  and  that  physical  modal  frequencies  and  damping  ratios  are  not  applicable.  A 
resulting  superstable  model  has  the  form  of  a  state-space  model  but  the  output  matrix  contains  the 
pulse  response,  i.e.,  the  output  matrix  contains  the  finite-time  system  dynamics.  The  system  ma¬ 
trix  contains  zeros  and  a  sub  diagonal  of  ones.  The  system  order  equals  the  number  of  time  steps. 
With  these  features  the  identification  is  direct,  the  system  is  stable,  and  simulations  have  the  fi¬ 
delity  of  the  pulse  response.  Here  we  show  that  a  superstable  model  derived  from  an  HPC- 
generated  wave-field  response  is  particularly  suited  to  large-domain  wave-propagation  simula¬ 
tions.  Our  example  is  an  FDTD  infrasonic  propagation  model. 

SUPERSTABLE  REPRESENTATION 

In  this  section  we  review  the  derivation  of  the  superstable  state-space  representation,8  and 
highlight  features  that  allow  direct  system  realization  from  finite-duration  large-wave-field  pulse 
responses.  State-space  equations  for  a  discrete  LTI  system  ( A,B,C )  are: 

x(k  + 1)  =  Ax(k)  +  Bu(k),  y{k)  =  Cx{k )  (2) 

where  x  is  an  n  -dimensional  state  vector;  u  is  an  r  -dimensional  input  vector;  and  y  is  a  q  - 
dimensional  output  vector.  A  is  the  n  x  n  system  matrix,  B  the  n  x  r  input  matrix,  and  C  the 
qxn  output  matrix,  k  is  the  time  index.  The  Markov  parameters  h,  ,  each  with  dimension  q x r  , 
characterize  the  pulse  response  of  this  system  as 

h{  =  CB,  h2  =  CAB,  h3  =  CA2B,...,  hk  =  CAk  lB,  k-l,2,...,p .  (3) 


The  realization  problem  is  stated  as  follows:  given  a  sequence  of  qxr  Markov  parameters,  find  a 
state-space  model  ( A,B,C )  such  that  Equation  (2)  holds  and  the  state-space  dimension  is  mini¬ 
mal.  49  A  minimum  realization  has  the  same  input-output  behavior  to  some  degree  of  accuracy.4 

Derivation 


To  derive  the  superstable  representation,  we  start  with  the  auto-regressive  moving-average 
model  with  exogenous  input  (i.e.,  ARX).  This  input-output  model  has  the  (order-p)  form 

y(k)  =  aly(k-l)  +  a2y(k-2)  +  ...+a  y(k- p)+ 

,  (4) 

[3\u(k  —  l)  +  (52u(k-2)  + ...+ (5pu(k- p). 


which  can  be  put  into  state-space  form  using 
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where  ( A,  BA')  is  the  observable  canonical  form. 
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With  the  Markov  parameters  in  Equation  (3),  the  unit  pulse  response  involving  the  p  Markov 
parameters,  starting  from  zero  initial  condition,  has  the  form 

y(k)  =  hlu(k-i)  +  h2u(k-2)  +  ...  +  hpu(k- p).  (6) 


Equation  (6)  is  the  finite  impulse  response  model.  If  we  view  Equation  (6)  as  a  special  case  of  the 
ARX  model  Equation  (4),  then  the  corresponding  version  of  Equation  (5)  is 
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Direct  substitution  reveals  that  ( AVBVC J  in  Equation  (7)  reproduces  the  p  Markov  parameters 
in  Equation  (3).  The  I  and  0  matrices  in  Equation  (7)  are  of  dimensions  qxq,  where  q  is  the 
number  of  outputs.  The  state-space  model  Equation  (7)  has  q  outputs,  r  inputs,  and  pq  states. 
For  a  system  with  a  small  number  of  inputs  and  a  large  number  of  outputs,  the  dimension  of  this 
model  is  large. 

Our  interest  is  a  minimal  system  for  the  large  output  problem,  and  so  we  will  work  with  the 
transposed  Markov  parameters  instead  of  the  original  Markov  parameters.  These  are 

hTk  =  (CAk~lB)T  =  BT(AT)k~{  CT  =  C*(A*)k~' B\  k  =  l,2,...,p.  (8) 


where  C  =  BT ,  A  =  A 7  ,  and  B  =Cr.  The  Markov  parameters  of  the  starred  system 
[a\b\c")  are  the  transposes  of  the  Markov  parameters  of  the  original  system  (A,  B,  c) .  Follow¬ 
ing  the  approach  above  leading  to  Equation  (7),  a  valid  state-space  realization  of  the  (a\b‘,C‘) 
system  is 
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{a;,b;,c;)  in  Equation  (9)  is  the  state-space  model  of  the  transposed  system  (a*,4,C*),  not  the 
original  system.  We  convert  (4, 4,4)  back  to  the  original  system.  In  so  doing,  we  arrive  at  the 
following  superstable  realization  for  the  original  system: 
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The  superstable  model  ( A2,B2,C2 )  has  q  outputs,  r  inputs,  and  p  states. 


Features  and  Significance 

A  superstable  model  is  a  finite-time  state-space  model  equivalent  to  the  finite  impulse  re¬ 
sponse  (FIR)  model.  Direct  substitution  reveals  that  Equation  (10)  reproduces  exactly  the  p 
Markov  parameters  in  Equation  (3).  In  applications  of  system  response  during  a  specific  duration 
of  interest,  this  representation  is  sufficient.  The  transient  dynamics  of  the  system  need  not  vanish 
during  this  period  for  the  model  to  be  applicable.  Although  the  representation  is  fundamentally 
finite-time,  the  model  may  still  be  used  to  predict  beyond  the  desired  duration  if  the  dynamics  of 
the  system  is  such  that  its  transient  response  is  negligible  at  the  end  of  the  finite  time  duration  of 
interest,  which  is  our  routine  practice. 

Interestingly  the  superstable  system  matrix  contains  no  dynamics.  The  system  matrix  has  zero 
eigenvalues  regardless  of  the  dynamics  of  the  system  that  it  models.  The  zero  eigenvalues  cause 
the  Markov  parameters  to  vanish  identically  beyond  a  finite  number  of  time  steps.  The  resulting 
stable  model  behavior  gives  the  technique  its  name.  In  a  typical  state-space  model,  the  system 
characteristic  information  (frequencies,  damping  ratios)  is  present  in  the  system  matrix,  which 
governs  the  transient  response  of  the  model.  The  concepts  of  modal  identification  are  not  applica¬ 
ble  to  wave  propagation  where  large  dimensions  preclude  standing  waves.  This  makes  the  super¬ 
stable  model  ideally  suited  to  external  wave  propagation,  since  while  the  waves  may  echo  and 
become  trapped  temporarily,  for  example  by  terrain  features  or  structures,  they  always  dissipate. 

As  explained  in  the  Introduction  we  select  our  FDTD  model  duration  to  have  negligible 
Markov  parameters  beyond  this  duration.  It  is  a  design  choice  to  assign  this  cutoff,  and  the  true 
FDTD  response  would  attenuate  toward  floating-point  precision.  This  is  analogous  to  a  decaying 
vibration  measurement.  Our  practice  is  to  impose  an  ending  taper  that  imposes  exactly  a  zero 
condition  at  the  end  of  the  FDTD  response,  and  consider  this  pulse  response  our  truth  model.  In 
the  case  of  a  source  with  excitation  beyond  the  duration  of  the  Markov  parameters,  the  model  will 
produce  a  response  with  accuracy  reflecting  the  design  dynamic  range.  However,  if  we  consider 
the  FDTD  computation  with  attenuation  toward  numerical  precision  to  be  the  truth  model,  then 
we  recognize  that  the  realization  is  exact  only  up  to  the  last  Markov  parameter  hp  and  not  beyond 

because  the  additional  Markov  parameters  formed  from  Equation  (10)  are  identically  zero.  This 
will  render  Equation  (10)  invalid  beyond  the  p  -time-step  duration,  which  highlights  the  impor¬ 
tance  of  selecting  the  duration  to  ensure  that  the  Markov  parameters  of  the  numerical  system  are 
negligible  from  that  point  on. 
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Figure  2.  FDTD  computation  of  superstable  model,  (a)  Areal  view  of  setting  (North  at  top)  showing 
source  location,  sensor  locations,  and  topography  and  partially  transparent  snapshot  overlays  from 
FDTD  output.  Snapshot  shows  waves  from  impulsive  source  passing  sensors  at  42.8  s  simulation 
time.  Insets  show  color  scales  for  topography  (lower  left)  in  meters,  and  pressure  (upper  right)  in 
Pascals,  (b)  Profiles  of  wave  speed  and  density,  (c)  Illustration  of  discrete  layers  used  to  approximate 
material  profiles,  (d)  Source  history,  (e)  Frequency  response  of  filter  used  to  generate  source. 
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In  our  wave  propagation  problems,  the  number  of  outputs  q  is  in  the  hundreds  of  thousands  or 
millions  whereas  the  number  of  inputs  r  is  1  for  a  single  source,  thus  a  huge  reduction  in  dimen¬ 
sions  is  achieved.  The  significance  is  that  a  time-consuming  realization  can  now  be  bypassed. 
Given  the  wave-field  Markov  parameters,  in  our  case  via  Fourier  analysis  to  calculate  the  pulse 
response,  a  minimum  superstable  model  of  the  system  can  be  immediately  obtained  without  any 
further  computation.  There  is  no  limitation  in  representing  any  propagation  or  other  physical  phe¬ 
nomenon  as  long  as  it  is  linear  and  the  time  duration  of  interest  produces  reasonable  dimensions 
for  computations,  although  state-space  based  model  reduction  techniques  can  be  applied  to  reduce 
the  dimension  of  the  model  further  if  warranted. 

INFRASONIC  PROPAGATION  MODEL 

Infrasound  is  sound  below  20  Hz.  It  can  propagate  very  long  distances  and  remain  at  measur¬ 
able  levels.  At  local  ranges,  e.g.,  on  the  order  of  10  kilometers,  terrain  relief  is  capable  of  scatter¬ 
ing  and  blocking  the  propagation.  This  is  because  the  infrasound  wavelengths  are  near  the  scales 
of  topographic  features.  These  interactions  have  been  observed  and  quantified  by  experimental 
data  and  models.10  Our  model,  which  we  created  for  FDTD  analyses  of  Equations  (1),  is  summa¬ 
rized  by  Figure  2.  It  spans  nearly  15  km  to  simulate  propagation  between  an  impulsive  just- 
above-ground  source  at  the  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC)  Big 
Black  Test  Site  (BBTS)  and  an  infrasound-sensing  array  at  the  ERDC  Waterways  Experiment 
Station  (WES).  Both  are  in  Mississippi.  Topographic  relief  between  these  locations  is  just  under 
70  m.  Figure  2  (a)  shows  the  areal  extent  of  the  model,  whose  grid  is  a  rectangular  box  in  three 
dimensions  with  a  height  of  3633  m.  The  grid  aligns  with  S  77.5°  E,  which  is  the  direction  be¬ 
tween  the  center  of  the  WES  array  labeled  “IML  1”  (IML  refers  to  the  locations  of  Inter- 
Mountain  Labs  sensors)  and  source  location  labeled  as  “BB  Source.”  Stretching  is  applied  so  that 
the  grid  has  higher  resolution  at  the  ground  surface  and  adjacent  to  the  source  and  sensor  loca¬ 
tions.  The  minimum  grid  spacing  is  2.75  m  horizontally  and  2.44  m  vertically.  The  grid  contains 
over  325  M  nodes. 

The  model  assigns  material  properties  at  the  appropriate  locations  to  represent  topographic  re¬ 
lief  and  atmospheric  profiles.  The  topography  is  a  stair-cased  interface  between  ground  and  air 
derived  from  digital  elevation  data.  Ground  properties  model  porous  sand;  these  are  flow  resis¬ 
tivity  a  =  5xl04  Pa-s/m2,  effective  density  p  =  11.7  kg/m3,  and  effective  wave  speed  c  =  151  m/s. 
The  model  uses  air-layer  properties  based  on  gridded  weather  forecast  data. 1  Our  interest  is  fore¬ 
cast  data  for  18:47  UTC  on  September  9,  2011,  which  is  the  time  of  a  10.9-kg  C-4  explosive  blast 
at  the  BBTS.  Figure  2  (b)  illustrates  the  atmospheric  conditions  and  wave  speed  and  density  pro¬ 
files  estimated  for  these  conditions.  The  model  includes  just  over  200  air  layers,  as  illustrated  in 
Figure  2  (c),  with  layer  thickness  increasing  with  elevation,  to  approximate  the  sound  speed  and 
density  profiles.  We  are  interested  in  numerical  accuracy  associated  with  10  nodes  per  wave¬ 
length.  The  fully  stretched  node  spacing  and  the  minimum  wave  speed  in  the  profile  results  in  a 
3 -Hz-accurate  bandwidth  for  the  grid.  This  is  high  enough  for  comparison  with  available  field 
data  from  the  blast  source. 


http://nationalmap.gov 
1  http://nomads.ncdc.noaa.gov 


Calibration  Analysis  to  Generate  Markov  Parameters 

The  first  FDTD  analysis  was  to  generate  response  data  for  calculating  Markov  parameters  of 
selected  output  domains.  The  analysis  applied  an  impulsive  dilatation-rate  source  at  coordinates 
within  the  model  corresponding  to  the  BBTS  blast  source  location,  1.2  m  above  the  ground  sur¬ 
face.  We  assigned  to  the  source  an  arbitrary  amplitude  (3160  m'/mVs  peak)  with  a  time  series  to 
ensure  that  the  frequency  content  of  the  source  energy  was  effectively  contained  within  the  3 -Hz 
model  bandwidth.  Figure  2  (d-e)  illustrate  the  source  time  series  and  frequency  content.  Our  the¬ 
ory  and  computations  are  for  linear  response,  and  thus  our  source  approximates  an  equivalent- 
linear  explosive  impulse. 

The  FDTD  analysis  ran  with  a  0.000225-s  sampling  interval  and  produced  512  output  time 
steps  with  a  0.1249-s  sampling  interval  (4.004-Hz  Nyquist  frequency),  reaching  a  level  with  ac¬ 
ceptable  decay  for  generating  Markov  parameters.  It  used  256  cores  for  12  h  45  min  on  a  Cray 
XE6.  The  cores  were  within  16  computational  nodes,  with  each  node  having  32  GB  of  memory. 
The  output  domain  presented  here  is  a  parallel-to-topography  slice  with  456xl544-output  loca¬ 
tions,  approximately  3  m  above  the  ground.  At  this  level  the  effect  of  terrain  relief  on  the  propa¬ 
gation  is  evident,  as  indicated  by  the  topography-imaging  character  of  Figure  3  and  subsequent 
overlay  results. 


Figure  3.  Areal  view  with  overlay  showing  pressure  level  from  FDTD  output  to  impulsive  source  of 
Figure  2  (d).  Level  is  for  sum-of-squares  of  signals  at  each  model  location.  Color  scale  (lower  left)  is 
in  dB  relative  to  1  Pa.  Second  overlay,  evident  around  the  source  location,  shows  pressure  level  of 
corresponding  superstable  model  in  dB  relative  to  the  FDTD  value — i.e.,  the  difference  in  dB  levels 
between  the  superstable  and  FDTD  integrated  signals.  This  overlay  is  set  to  be  transparent  below 
±0.005  dB,  an  arbitrary  threshold  that  is  approximately  0.06  percent  relative  error.  It  saturates  the 

color  scale  at  +0.05  dB. 

Figure  3  presents  a  measure  of  pressure  level  throughout  the  aboveground  domain  from  the 
FDTD  simulation.  From  the  pulse  response  data  underlying  the  image  in  Figure  3,  we  created  a 
512-state  superstable  model  by  populating  the  C  Matrix  of  Equation  (10).  As  expected,  the  fidel- 


http://www.erdc.hpc.mil/hardware/index.html 
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ity  of  the  superstable  model  to  the  FDTD  simulation  is  high  when  excited  using  the  calibration 
source,  which  was  down  sampled  to  the  FDTD-output  sampling  interval  (i.e.,  the  discrete  signal 
outlined  by  the  asterisks  in  Figure  2  (d)).  The  difference  in  superstable  and  FDTD  pressure  levels, 
also  illustrated  in  Figure  3,  exceeds  0.005  dB  only  in  the  vicinity  of  the  source. 


Figure  4.  FDTD  and  superstable  model  signals  from  3  m  above  the  ground  compared  with  measured 
ground-level  data  from  experimental  blast.  The  lower  two  signals  are  from  locations  corresponding 
to  the  nearby  IML  1  and  2  sensors,  as  indicated  in  Figure  2,  and  the  upper  signals  are  from  the  loca¬ 
tion  of  the  IML  sensor  to  the  South.  The  blue  lines  are  the  experimental  data  filtered  to  a  3-Hz 
bandwidth.  These  lines  show  ambient  noise  as  well  as  arriving  energy  from  the  blast.  The  red  lines 
are  from  the  superstable  model.  Green  lines,  which  are  hidden  exactly  below  the  red  lines,  are  FDTD 
signals.  Peak-to-peak  values  are  listed  for  each  location  in  order  from  IML  1.  The  inset  shows  a  cross 
correlation  result  between  the  source  and  FDTD  IML-1  signals,  indicating  that  the  arrival  is  close  to 
the  reported  experimental  arrival,  which  was  42.7  s. 

With  the  calibration-model  energy  levels  as  background.  Figure  4  graphs  field  signals  from 
the  IML  locations  with  corresponding  model  signals.  The  IML  sensors  have  a  30-Hz  nominal 
range.  We  filtered  the  measured  signals  with  a  3-Hz  cutoff  and  down  sampled  to  a  0.1240-s  sam¬ 
pling  rate  to  allow  comparison.  The  model  signals  show  a  waveform  similar  to  the  field  signals  at 
the  time  of  arrival  of  the  source  impulse.  The  peak-to-peak  amplitude  values  given  in  Figure  4 
suggest  that  a  peak  dilatation  rate  close  to  20  nr/nr/s  for  the  applied  source  signal  in  Figure  2  (d) 
would  have  produced  a  good  amplitude  match.  The  experimental  data  shows  considerable  noise, 
relative  to  the  source  signal,  and  the  source  arrival  appears  to  drop  below  noise  levels  within 
about  2.5  s.  The  model  data  matches  this  decay  well.  The  suggestion  is  that  the  3-Hz-bandwidth 
impulsive-source  FDTD  model  has  good  fidelity  with  the  experimental  data,  and  thus  so  do  the 
Markov  parameters. 
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Figure  5.  Results  from  a  sequence  of  filtered  source  pulses,  generated  from  superimposed  2/3-,  4/3-, 
and  6/3-Hz  sinusoids  over  a  finite  duration,  (a)  Areal  view  of  FDTD-output  snapshot  at  42.8  s.  Inset 
shows  graph  of  source  history,  (b)  Oblique  view,  looking  toward  the  Northwest,  with  overlay  of  pres¬ 
sure  level  from  superstable  model  using  down-sampled  FDTD  source.  Level  is  for  sum-of-squares  of 
signals  at  each  model  location.  Color  scale  (left)  is  in  dB  relative  to  1  Pa.  Second  overlay  shows  pres¬ 
sure  level  of  the  superstable  model  relative  to  the  FDTD  level.  It  is  transparent  where  below  a  ±0.005 

dB  threshold. 
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Figure  6.  Effect  of  model-order  reduction  on  results  using  sequence  of  finite-duration  source  pulses, 
(a)  Singular  value  plots  from  each  of  192  data  blocks  in  output,  (b)  512-state  superstable  model.  View 
toward  the  Southeast  shows  overlay  of  pressure  level  in  dB  relative  to  1  Pa  calculated  using  sum-of- 
squares  of  signals  at  each  model  location.  Color  scale  is  to  the  left.  Second  overlay  shows  pressure 
level  of  the  superstable  model  relative  to  the  FDTD  level.  Color  scale  is  to  the  right.  Here  the  trans¬ 
parent  threshold  is  ±0.05  dB,  or  approximately  0.6%  relative  error,  (c-d)  384-  and  320-state  reduced- 
order  models,  respectively,  (e)  Signals  at  the  IML  locations  from  FDTD  model  (blue),  512-state  su¬ 
perstable  model  (green),  384-state  reduced  model  (red),  and  320-state  reduced  model  (turquoise). 
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Verification  and  Performance  of  Superstable  Model 

The  second  FDTD  analysis  was  to  verify  that  a  simulation  using  the  superstable  model  was 
accurate  for  an  independent  source.  The  analysis  applied  a  sequence  of  source  pulses,  generated 
from  superimposed  10-s-duration  2/3-,  4/3-,  and  6/3-Hz  equal-amplitude  sinusoids,  at  the  location 
of  the  source  in  the  calibration  model.  We  filtered  the  pulses  with  the  response  in  Figure  2  (e)  to 
ensure  the  source  energy  was  within  the  bandwidth  of  the  superstable  model.  The  peak  amplitude 
was  nominally  20  m3/m3/s.  Figure  5  illustrates  the  source  time  series,  a  snapshot  of  the  propaga¬ 
tion,  and  the  wave-field  pressure-level  response  superimposed  on  the  topography  with  elevation 
exaggerated  3x  for  viewing.  This  FDTD  analysis  also  processed  512  output  time  steps  with  a 
0.1249-s  sampling  interval  and  used  nearly  the  same  computational  duration  on  256  cores  of  the 
Cray  XE6  as  the  calibration  analysis,  i.e.,  12  h  45  min. 

Figure  5  (b)  illustrates  the  fidelity  of  a  corresponding  512-step  superstable  simulation  using 
this  second  source.  The  pressure  level  relative  to  the  FDTD  simulation  exceeds  0.005  dB  in  rings 
surrounding  the  source,  and  0.05  dB  in  the  immediate  vicinity  of  the  source,  which  are  very  small 
differences  in  a  practical  context.  The  superstable  state-update  and  output  time-series  calculations 
with  save-to-file  operations  lasted  170  s.  These  used  MATLAB®-coded  software  on  a  2.66  GHz 
Apple  MacBook  Pro  5,3.  The  software  accessed  data  using  memory-mapping  capabilities  to 
avoid  reading  the  Equation-(IO)  C  -matrix  data  into  memory.  This  data  resided  in  a  single¬ 
precision  1.44-GB  file  on  an  external  drive  connected  by  an  800  Mb/sec  serial  interface.  The  C  - 
matrix  comprised  512-step  time  series  of  704064  outputs.  Considering  number  of  cores,  average 
memory  per  core,  and  analysis  duration,  the  superstable  simulation  used  approximately  one  mil¬ 
lionth  of  the  computation  resources  of  the  FDTD  model.  While  this  comparison  is  for  an  output 
that  is  only  about  0.2  percent  of  the  FDTD  computational  nodes,  the  savings  clearly  shows  the 
efficiency  in  reuse  of  pre-selected  wave-field  results  by  designing  the  HPC  analysis  as  a  supersta¬ 
ble-model  pre-processing  step. 

Reduced-Order  Models  vs.  Numerical  Prediction 

Our  superstable  state-space  dynamic  calculations  are  a  streamlined  version  of  the  Equation-(2) 
update  equations  by  virtue  of  the  simple  forms  of  matrices  A  and  B  in  Equation  10.  The  states 
update  from  the  source  amplitudes  in  a  loop  over  the  time  steps,  and  the  output  is  the  Markov  pa¬ 
rameter  time  series  multiplied  by  the  state  time  series.  Using  this  calculation  approach  is  a  matter 
of  preference.  In  state-space  form  the  model  is  compatible  with  available  state-space  based  com¬ 
putational  and  analysis  tools  including  model  reduction,  model  inverse,  and  Kalman  filtering. 
These  tools  are  not  always  available  for  the  unit  pulse  response  model  that  the  superstable  model 
replaces.  Here  we  apply  model  reduction  to  the  superstable  model  and  illustrate  the  error  that  re¬ 
sults. 

Again  using  the  second  FDTD  response  to  the  Figure  5  source  time  series  as  the  standard,  the 
512-state  superstable  model  was  reduced  using  MATFAB®  algorithms  for  balanced  realizations 
and  model  reductions,  and  used  in  simulations  performed  with  the  corresponding  down  sampled 
source.  Because  of  memory  restrictions  on  the  MacBook  Pro,  we  performed  the  balancing  and 
model  reductions  in  a  series  of  192  blocks  of  data,  each  with  3667  outputs.  (In  trial  analysis  we 
used  as  few  as  8  blocks  of  88008  outputs,  but  arbitrarily  chose  the  larger  number  of  blocks  in 
these  calculations).  Figure  6  (a)  shows  the  singular-value  plots  for  these  blocks,  with  the  differ¬ 
ence  between  ends  of  the  model  indicated.  Based  on  inspection  of  this  graph,  we  produced  mod¬ 
els  with  384  and  320  states  to  examine  the  increasing  error. 

Figure  6  (b-d)  show  the  overlays  of  pressure  level  from  the  512-,  384-,  and  320-state  models, 
respectively.  The  view  is  toward  the  Southeast  to  highlight  the  errors  in  the  reduced-order  models 
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at  the  far  end  of  the  model  from  the  source.  Here  the  difference  in  pressure  levels  between  the 
512-state  superstable  model  and  the  FDTD  result  does  not  rise  above  a  ±0.05  dB  threshold,  or  0.6 
percent  relative  error.  In  contrast,  both  the  384-  and  320-state  models  show  errors  that  pass  this 
threshold,  with  the  latter  encompassing  considerable  area. 

This  relative  error,  however,  is  insignificant  when  considering  the  FDTD  propagation  system 
and  the  level  of  fidelity  it  has  to  real  signals,  for  example  those  in  Figure  4.  This  is  underscored 
by  the  signals  comparison  in  Figure  6  (e),  which  shows  FDTD  model  signals  and  the  512-,  384-, 
and  320-order  state-space  model  signals  from  the  locations  of  the  IML  sensors.  A  difference  is 
just  noticeable  with  the  line  thickness  used  in  the  graph  and  the  significant  digits  used  in  the 
peak-to-peak  amplitude  values.  (The  signals  plot  in  the  order  listed  by  the  Figure  6  (e)  caption, 
and  are  very  close  to  one  another,  so  that  the  turquoise  and  red  lines  are  all  that  are  visible).  The 
implication  is  that  model-order  reduction  can  serve  a  useful  purpose  when  stored  file  size  is  of 
importance  and  the  reduced  accuracy  is  acceptable.  In  the  384-  and  320-state  reduced-order  cases, 
the  stored  systems  would  benefit  from  a  reduced  size  C  -matrix,  75  percent  and  62  percent  of  the 
original  1.44-GB  size,  respectively,  while  requiring  slightly  more  storage  for  the  non-trivial  A 
and  B  matrices. 

CONCLUSION 

This  paper  reviews  and  demonstrates  the  discrete-time  superstable  model.  The  model  was  mo¬ 
tivated  by  application  of  the  Eigensystem  Realization  Algorithm  to  wave-field  systems  from  HPC 
analyses.  Applications  of  ERA  revealed  system  features  allowing  superstable  model  assignment. 
Most  importantly  the  pulse  response  and  its  decay  over  the  domain  are  captured  in  a  finite  dura¬ 
tion,  and  decay  to  zero  beyond  a  finite  number  of  time  steps  implies  a  system  with  zero  eigenval¬ 
ues.  This  makes  the  superstable  model  ideally  suited  to  wave  propagation  where  large  dimensions 
preclude  standing  waves  and  wave  energy  always  dissipates.  The  model  is  most  efficient  when 
the  number  of  input  sources  is  small,  the  number  of  outputs  is  very  large,  and  the  number  of  sam¬ 
ples  associated  with  the  finite-time  duration  of  interest  is  small.  When  Markov  parameters  are 
computed  from  HPC-derived  input-output  data,  a  superstable  model  can  be  assigned  immediately 
without  additional  computation.  While  a  superstable  model  is  a  finite-time  state-space  model 
equivalent  to  the  finite  impulse  response  model,  in  state-space  form  the  model  is  compatible  with 
a  wide  variety  of  advantageous  state-space  based  computational  and  analysis  tools.  These  tools 
include  dynamic  simulation,  model  reduction,  model  inverse,  and  Kalman  filtering. 

Our  work  here  demonstrated  propagation-system  identification  with  pulse  response  data  from 
an  infrasonic  propagation  analysis  in  which  terrain  affects  the  character  of  the  pressure  level 
throughout  the  wave  field.  We  showed  the  high  fidelity  of  a  superstable  simulation  to  a  field- 
validated  supercomputer  analysis  of  propagation  from  an  impulsive  source,  and  verified  its  accu¬ 
racy  with  a  numerical  analysis  using  a  different  source.  The  superstable  result,  performed  on  a 
laptop,  used  approximately  one  millionth  of  the  computation  resources  of  the  supercomputer 
analysis,  showing  the  benefit  of  designing  an  HPC  analysis  as  a  superstable-model  pre-processing 
step.  Application  of  balanced  realization  and  model-order  reduction  algorithms  to  the  superstable 
model  demonstrated  the  advantage  of  preparing  a  model  in  state-space  form  when  accuracy  is 
acceptable  and  reduced-size  file  storage  is  required.  We  conclude  that,  by  superstable-model 
identification,  we  are  able  to  create  a  reusable  and  reducible  state-space  propagation-system 
model  for  a  given  source  configuration  that  accurately  simulates  large-wave-field  time-domain 
propagation  using  a  fraction  of  the  original  computational  resources. 
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